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Human learning is a complex phenomenon requiring flexibility to 
adapt existing brain function and precision in selecting new neuro- 
physiological activities to drive desired behavior [1, 2]. These two 
attributes - flexibility and selection - must operate over multiple 
temporal scales as performance of a skill changes from being slow 
and challenging to being fast and automatic [3, 4]. Such selec- 
tive adaptability is naturally provided by modular structure [5, 6, 7], 
which plays a critical role in evolution [8], development [6], and 
optimal network function [7, 9]. Using functional connectivity mea- 
surements of brain activity acquired from initial training through 
mastery of a simple motor skill, we investigate the role of modu- 
larity in human learning by identifying dynamic changes of modular 
organization spanning multiple temporal scales. Our results indi- 
cate that flexibility, which we measure by the allegiance of nodes to 
modules, in one experimental session predicts the relative amount 
of learning in a future session. We also develop a general statistical 
framework for the identification of modular architectures in evolv- 
ing systems, which is broadly applicable to disciplines where network 
adaptability is crucial to the understanding of system performance. 
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The brain is a complex system, composed of many in- 
teracting parts, which dynamically adapts to a continually 
changing environment over multiple temporal scales. Over 
relatively short temporal scales, rapid adaptation and contin- 
uous evolution of those interactions or connections form the 
neurophysiological basis for behavioral adaptation or learning. 
At small spatial scales, stable neurophysiological signatures of 
learning have been best demonstrated in animal systems at 
the level of individual synapses between neurons [10, 11, 12]. 
At a larger spatial scale, it is also well known that specific 
regional changes in brain activity and effective connectivity 
accompany many forms of learning in humans — including the 
acquisition of motor skills [1, 2]. 

Learning-associated adaptability is thought to stem from 
the principle of cortical modularity [13]. Modular, or nearly 
decomposable [5], structures are aggregates of small subsys- 
tems (modules) that can perform specific functions without 
perturbing the remainder of the system. Such structure pro- 
vides a combination of compartmentalization and redundancy, 
which reduces the interdependence of components, enhances 
robustness, and facilitates behavioral adaptation [6, 14]. Mod- 
ular organization also confers evolvability on a system by re- 
ducing constraints on change [6, 7, 15, 16]. Indeed, a putative 
relationship between modularity and adaptability in the con- 
text of human neuroscience has recently been posited [17, 18]. 
To date, however, the existence of modularity in large-scale 
cortical connectivity during learning has not been tested di- 
rectly. 

Based on these aforementioned theoretical and empirical 
grounds, we hypothesized that the principle of modularity 
would characterize the fundamental organization of human 



brain functional connectivity during learning. More specifi- 
cally, based on several studies relating the neural basis of mod- 
ularity to the development of skilled movements [19, 20, 21], 
we expected that functional brain networks derived from ac- 
quisition of a simple motor skill would display modular struc- 
ture over the variety of temporal scales associated with learn- 
ing [4]. We also hypothesized that that modular structure 
would change dynamically during learning [1, 3], and that 
characteristics of such dynamics would be associated with 
learning success. 

We tested these predictions using functional magnetic res- 
onance imaging (fMRI), an indirect measure of local neuronal 
activity [22], in healthy adult subjects during the acquisition 
of a simple motor learning skill composed of visually cued fin- 
ger sequences. We derived low-frequency (0.06-0.12 Hz) func- 
tional networks from the fMRI data by computing the tempo- 
ral correlation between activity in each pair of brain regions to 
construct weighted graphs or whole-brain functional networks 
[23, 24, 25] (See Supplementary Information for details). This 
network framework enabled us to estimate a mathematical 
representation of modular or community organization, known 
as 'network modularity', for each individual over a range of 
temporal scales. We evaluated the evolution of network con- 
nectivity over time using a novel mathematical framework 
[28], and we tested its relationship with learning using the 
Pearson correlation coefficient. See Methods for details of the 
sample, experimental paradigm, and methods of analysis. 



Results 

Static Modular Structure. We investigated network organi- 
zation over multiple temporal scales — over days, hours, and 
minutes — during motor learning [3, 4] (see Figure IB), and we 
found it to be significantly modular at each scale (see Figure 
2A-C). A diagnostic measure of the amount of network mod- 
ularity in the system is the modularity index Q (See Methods 
for a mathematical definition). Importantly, while an increase 
in the modularity index Q indicates a significant segregation 
of system structure into distinct modules or communities, the 
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number of modules that we uncovered is significantly smaller 
than that expected in a random system, indicating a comple- 
mentary integration. The scaling of architectural organization 
over multiple temporal resolutions is consistent with previous 
findings for the human cortex in both the frequency [26] and 
spatial domains [9, 27]. Furthermore, the temporal structure 
of this architecture is graded in the sense that fewer modules 
(about 3) on longer time scales (Figure 2A-B) are comple- 
mented by more modules (about 4) on shorter time scales 
(Figure 2C). This graded structure is analogous to that found 
in the nested modular networks of underlying brain anatomy 
where few modules uncovered at large spatial scales are com- 
plemented by more modules at smaller spatial scales [9]. 

Dynamic Modular Structure. We next consider evolvability, 
which is most readily detected when the organism is under 
stress [8] or when acquiring new capacities such as during 
external training in our experiment. We found that the com- 
munity organization of brain connectivity reconfigured adap- 
tively over time. Using a recently-developed mathematical 
formalism to assess the presence of dynamic network recon- 
figurations [28], we constructed multilayer networks in which 
we link the network for each time window (see Figure 3A) 
to the network in the time windows before and after (Figure 
3B) by connecting each node to itself in the neighboring win- 
dows. We then measured modular organization [53, 54, 29] on 
this linked multilayered network in order to find long-lasting 
modules [28]. 

To verify the reliability of our measurements of dynamic 
modular architecture, we introduced three null models based 
on permutation testing (Figure 3C). We found that cortical 
connectivity is specifically patterned, which we concluded by 
comparison to a 'connectional' null model in which we scram- 
bled links between nodes in each time window [30]. Further- 
more, cortical regions maintain these individual connectivity 
signatures that define community organization, which we de- 
termined by comparison to a 'nodal' null model in which we 
linked a node in one time window to a randomly chosen node 
in the previous and next time windows. Finally, we found 
that functional communities display a smooth temporal evo- 
lution, which we identified by comparing diagnostics derived 
from the true multilayer network structure to those derived 
from a temporally permuted version (see Figure 3D). We con- 
structed this 'temporal' null model by randomly re-ordering 
the multilayer network layers in time. 

By comparing the structure of the cortical network to that 
of the null models, we found that the human brain displayed 
a heightened modular structure in which more modules of 
smaller size were discriminable as a consequence of the emer- 
gence and extinction of modules in cortical network evolution. 
The 'stationarity' of communities, which is defined as the aver- 
age correlation between partitions over consecutive time steps 
[31], was also higher in the human brain than in the con- 
nectional or nodal null models, indicating a smooth temporal 
evolution. 

Learning. Given the dynamic architecture of brain connectiv- 
ity, it is interesting to ask whether the specific architecture 
changes with learning — either at a gross scale through an 
adaptation in the number or sizes of modules or at a finer 
scale through alterations in the nodal composition of mod- 
ules. Empirically, we found no significant differences between 
experimental sessions in the coarse diagnostics. To quantify 
finer-scale architectural fluctuations, we introduced the no- 
tion of node flexibility, which we define using the multilayer 
framework to be the number of times that each node changes 



module allegiance, normalized by the total possible number of 
changes (see Supplementary Materials for details). The flex- 
ibility of the network as a whole is then defined as the mean 
flexibility over all nodes. 

Network flexibility is a biologically meaningful measure 
that captures changes in the local properties of individual 
network elements. We found that network flexibility changed 
during the learning process - first increasing and then decreas- 
ing (see Figure 4A). In particular, the flexibility of a partici- 
pant in one session could be used as a predictor of the amount 
of learning (as measured by the improvement in the time re- 
quired to complete the sequence of motor responses) in the 
following session (Figure 4B). Regions of the brain that were 
most responsible for this predictive power of individual differ- 
ences in learning were distributed throughout the cortex, with 
strong loadings in the medial, inferior, and pre-frontal, presup- 
plementary motor, posterior parietal, and occipital cortices 
(Figure 4C-D). The ability to predict future learning capacity 
could not be reliably derived from fMRI activation, support- 
ing our conclusion that flexibility provides a useful approach 
for modeling system evolvability. 

Our results indicate that flexibility is sensitive to both 
intra-individual and inter-individual variability. Across par- 
ticipant, we found that network flexibility was modulated by 
learning (Figure 4A). However, we also found that each par- 
ticipant displayed a characteristic flexibility. The variation in 
flexibility over participants was larger than the variation in 
flexibility across sessions, as measured by the intra-class cor- 
relation coefficient (ICC « 0.56; F-statistic F(17,34) « 4.85, 
p^4 x 10~ 5 ). 



Discussion 

Modularity of Functional Connectivity. Modularity is an in- 
tuitively important property for dynamic, adaptable systems. 
The accompanying system decomposability provides the nec- 
essary structure for complex reconfigurations. Modularity can 
be a property of morphology, as has been widely described in 
the context of evolution and development [8, 15, 16], as well 
as of the interconnection patterns of social, biological, and 
technological systems [53, 54]. More pertinent to this pa- 
per, recent evidence suggests that modular organization over 
several spatial scales, or hierarchical modularity, also charac- 
terizes the large-scale anatomical connectivity of the human 
brain [9, 27], as well as the spontaneous fluctuations [40, 41] 
thought to stem from anatomical patterns [42]. However, the 
putative relationship between adaptability and modular struc- 
ture has not been previously explored in the context of the 
brain connectome. 

In the present study, we have shown that the functional 
connectivity of the human brain during a simple learning 
paradigm is non- homogenous. Rather, it is segregated into 
communities that can each perform unique functions. This 
segregation of connectivity structure was expressed consis- 
tently over the scale of days, hours, and minutes, suggesting 
that community structure provides a generalizable framework 
for the successful actuation of temporally distinct phenomena 
[16]. However, it is also notable that connectivity at the short- 
est temporal scale displayed higher variability, perhaps re- 
flecting the necessity for dynamic modulation of human brain 
function over relatively short intervals during learning [3] . In 
light of historically strict definitions of cognitive modules as 
completely encapsulated structures [43] , it is important to em- 
phasize that the modules that we have uncovered here remain 
integrated with one another by a complex pattern of weak 
interconnections. 
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Dynamic Network Evolution. Efforts to characterize both rest- 
ing state [44] and task-based large-scale connectivity of human 
brain structure and function [23, 24, 25] have focused almost 
exclusively on static representations of underlying connectiv- 
ity patterns. However, both scientific intuition and recent evi- 
dence suggest that connectivity can be modulated both spon- 
taneously [32] and by exogenous stimulation [1]. The explo- 
ration of temporally evolving network architecture therefore 
forms a critical new frontier in neuroscience. 

Our exploration of dynamic community structure in an ex- 
perimental paradigm that requires neurophysiological adapt- 
ability provides insight into the organizational principles sup- 
porting successful brain dynamics. Similar to social systems 
[31], we found that community organization changed smoothly 
with time, displaying coherent temporal dependence on what 
had gone before and what came after, a characteristic com- 
patible with complex long-memory dynamical systems [45]. 

In addition to global adaptability, we found that diverse 
regions of the brain performed different roles within those 
communities: Some maintain community allegiance through- 
out the experiment (low flexibility nodes), and others con- 
stantly shift allegiance (high flexibility nodes). Biologically, 
this network flexibility might be driven by physiological pro- 
cesses that facilitate the participation of cortical regions in 
multiple functional communities. Learning a motor skill in- 
duces changes in both the structure and connectivity of the 
cortex[33, 34], which is accompanied by increased excitability 
and decreased inhibition of neural circuitry [35, 36, 37]. How- 
ever, it is plausible that flexibility might also be driven by 
task-dependent processes that require the capacity to balance 
learning across subtasks. For example, the particular experi- 
ment utilized in this study demanded that subjects master the 
performance of using the response box, decoding the stimu- 
lus, performing precise movements, shifting attention between 
stimuli, and switching between different movement patterns. 



Flexibility and Learning. Importantly, the inherent temporal 
variability in network structure measured by nodal flexibility 
was not a stable signature of an individual's functional orga- 
nization but was instead modulated by consecutive stages of 
learning - first increasing and then decreasing as movement 
time stabilized in the later stages of learning [3]. The modula- 
tion of flexibility by learning was evident not only at the group 
level but also in individuals. The amount of flexibility in each 
participant could be used to predict that participant's learn- 
ing in a following experimental session. In addition to sup- 
porting the theoretical utility of accessible but often ignored 
higher-order (bivariate, multivariate) statistics of brain func- 
tion, this result could potentially be used to inform decisions 
on how and when to train individuals on new tasks depending 
on the current flexibility of their brain. From this work alone, 
however, we are unable to determine whether or not learning 
is the only possible modulator of flexibility. Complementary 
experiments could be designed to test whether flexibility is 
also modulated by fatigue or exogenous stimulants in order 
to increase subsequent skilled learning. We also found that 
inter-individual variability in flexibility was larger than intra- 
individual variability, indicating that flexibility might be a re- 
liable indicator of a given subject's brain state. Consequently, 
our methodology could potentially be of use in predicting a 
given individual's response to training or neurorehabilitation 
[38, 39]. 

Flexibility might be a network signature of a complex un- 
derlying cortical system characterized by noise [46]. Such a 
hypothesis is bolstered by recent complementary evidence sug- 
gesting that variability in brain signals also supports mental 
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effort in a variety of cognitive operations [47], presumably by 
aiding the brain in switching between different network con- 
figurations as it masters a new task. Indeed, the theoretical 
utility of noise in a nonlinear dynamical system like the brain 
[48] lies in its facilitation of transitions between network states 
or system functions [49] and therefore helps to delineate the 
system's dynamic repertoire [50]. However, despite the ap- 
parent plausibility that network flexibility and cortical noise 
are related, future studies are necessary to directly test this 
hypothesis. 

Methodological Considerations. The construction of brain 
networks from continuous association matrices, such as those 
based on pairwise correlation or coherence, has historically 
been performed by applying a threshold to the data to con- 
struct a binary graph in which an edge exists if the association 
between the nodes it connects was above the threshold and 
does not exist otherwise [23, 24, 25]. However, the statistical 
validity of that method is hampered by the need to choose 
an arbitrary threshold as well as by the discretization of in- 
herently continuous edge weights. In the current work, we 
have instead used fully weighted networks in which connec- 
tions retain their original association value unless that value 
was found to be insignificant (based on statistical testing em- 
ploying a false discovery rate (FDR) correction for multiple 
comparisons) [52]. Future studies comparing multiple network 
construction techniques will be important to statistically as- 
sess the added value of weighted edge retention in the assess- 
ment of network correlates of cognition. 

Second, partitioning a set of nodes into a set of communi- 
ties is NP-hard [55] , meaning that modularity algorithms pro- 
duce many near-optimal partitions of the network [56]. The 
number of near-optimal partitions is higher in large networks, 
or those with edges weighted with only l's and 0's [56]. Here 
we have chosen to study small weighted networks in which 
the number of near-optimal partitions is small. Nevertheless, 
we have systematically explored the partition landscape by 
iterative optimization. Accordingly, we report mean modu- 
larity estimates that our results suggest are representative of 
the system (see Supplementary Materials for further details). 
However, further work is necessary to measure common com- 
munity assignments within the ensemble of partitions in or- 
der to identify consistently segregated groups of brain regions. 
This will aid in further exploration of the biological relevance 
of the uncovered communities. 

Finally, the statistical validation of community structure 
in social and biological systems is complicated by several fac- 
tors. For example, many investigations, especially in social 
systems, are hindered by their small number of instantiations. 
In our work, the relatively large number of subjects in con- 
junction with estimations of multiple networks over various 
temporal scales facilitated a stringent statistical assessment 
of community structure both in comparison to randomly con- 
nected graphs and, as we have developed for dynamic net- 
works, to graphs where nodal identities or times were scram- 
bled. An important future area of research will focus on the 
development of alternative null models that are not perfectly 
random but which assume increasingly biologically realistic 
network architectures. 



Conclusion 

Consistent with our hypotheses, we have identified signifi- 
cant modular structure in human brain function during learn- 
ing over a range of temporal scales: over days, hours, and 
minutes. Modular organization over short temporal scales 
changed smoothly, suggesting system adaptability. The com- 
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position of functional modules displayed temporal flexibility 
that was modulated by early learning, varied over individuals, 
and was a significant predictor of learning in subsequent ex- 
perimental sessions. Furthermore, we developed and reported 
a general framework for statistical validation of dynamic mod- 
ular architectures in arbitrary systems. More generally, our 
evidence for adaptive modular organization in global brain ac- 
tivity during learning provides critical insight regarding the 
dependence of system performance on underlying architec- 
ture. 

Materials and Methods 

Twenty-five right-handed participants (16 female, 9 male; mean age 24.25 years) vol- 
unteered with informed consent in accordance with the UCSB internal Review Board. 
After exclusions for task accuracy, incomplete scans, and abnormal Magnetic Res- 
onance Imaging (MRI), 18 participants were retained for subsequent analysis. All 
participants had less than 4 years of experience with any one musical instrument, 
had normal vision, and had no history of neurological disease or psychiatric disorders. 
Participants were paid for their participation. 

The experimental framework consisted of a simple motor learning task in which 
subjects responded to a visually cued sequence by generating responses using the 
4 fingers of their non-dominant hand (thumb excluded) on a custom response box. 
Participants were instructed to respond swiftly and accurately. Visual cues were pre- 
sented as a series of musical notes on a 4-line music staff such that the top line of 
the staff mapped to the leftmost key depressed with the pinkie finger. Each 12-note 
sequence contained 3 notes per line, which were randomly ordered without repetition 
and free of regularities such as trills and runs. The number and order of sequence 
trials was identical for all participants. All participants completed 3 training sessions 
in a 5-day period, and each session was performed inside the MRI scanner. 

Functional MRI (fMRI) recordings were conducted using a 3.0 T Siemens Trio 
with a 12-channel phased-array head coil. For each functional run, a single-shot echo 
planar imaging that is sensitive to blood oxygen level dependent (BOLD) contrast 
was used to acquire 33 slices (3 mm thickness) per repetition time (TR), with a TR 
of 2000 ms, an echo time of 30 ms, a flip angle of 90 degrees, a field of view of 192 
mm, and a 64 X 64 acquisition matrix. Image preprocessing was performed using 
the FMRIB (Oxford Centre for Functional Magnetic Resonance Imaging of the Brain) 
Software Library (FSL), and motion correction was performed in MCFLIRT (Motion 
Correction using FMRIB's Linear Image Registration Tool). Images were high-pass 
filtered with a 50 s cutoff period. Spatial smoothing was performed using a kernel 
where full width at half maximum was 8 mm. Signals were normalized globally to 
account for transient fluctuations in intensity. 

The whole brain is parcellated into a set of N regions of interest that corre- 
spond to the 112 cortical and subcortical structures anatomically identified in FSL's 
Harvard-Oxford atlas. For each individual fMRI data set, we estimate regional mean 
BOLD time series by averaging voxel time series in each of the N regions. These 
regional time series are then subjected to a wavelet decomposition to reconstruct 
wavelet coefficients in the 0.06-0.12 Hz range (scale two). We estimate the correla- 
tion or coherence Aij between the activity of all possible pairs of regions % and j to 
construct N X TV functional connectivity matrices A (see Figure 1A). Individual 
elements of Aij are subjected to statistical testing, and the value of all elements 
that do not pass the false discovery rate correction for multiple comparisons are set 
to zero; otherwise, the values remain unchanged. The complete set of weighted net- 
work nodes is partitioned into communities by maximizing the modularity index Q 
with respect to the connectivity of a random null model[53, 54]. In the simplest static 
case, supposing that node % is assigned to community Qi and node j is assigned to 



community Qj, then Q is defined as 

Q = 52[Ai3-Pij]8(gi,9i), m 

ij 

where $(Qi->Qj} — 1 if = Qj and it equals otherwise, and Pij is the 
expected weight of the edge connecting node I and node J under a specified null 
model. (Note: a more complex formula is used in the dynamic network case; see 
Supplementary Information.) The elements of the matrix Aij are weighted by the 
functional association between regions and we thoroughly sample the distribution of 
partitions that provide near-optimal Q values [56]. The functional connectivity is 
termed 'modular' if the value of Q is larger than that expected from random network 
null models that control for both the mean and variability of connectivity. 

We tested for static modular structure on these individual networks and on 
dynamic network structure on a multi-network object created by linking networks 
between time steps [28]. In both cases, we assess modular organization using the 
modularity Q and the number of modules U. In the dynamic case, we also used two 
additional diagnostics to characterize modular structure including the mean module 
size S and the stationarity of modules £. We defined S to be the mean number of 
nodes per community over all time windows over which the community exists. We 
used the definition of module stationarity from Ref. [31]. We started by calculating 
the autocorrelation function U(t) of two states of the same community G(t) at 
t = 1 time steps apart using the formula 



U(t) 



\G(t )nG(t + t)\ 
\G(t )uG(t + t)\ 



[2] 



where to ' s the time at which the community is born, If? (to) H G(to + t)| 
is the number of nodes that are members of both G(to) ar| d G(to + t), and 
\G(t ) U G(to + t)| is the total number of nodes in G(tn) U G(t + t) [31]. 
We defined t' to be the final time step before the community is extinguished. The 
stationarity of a community is then 



Ettl u(t,t + i) 



t'-to-l 

which is the mean autocorrelation over consecutive time steps [31]. 



[3] 



In principle, modular architecture might vary with learning by displaying changes 
in global diagnostics such as the number of modules or the modularity index Q or by 
displaying more specific changes in the composition of modules. To measure changes 
in the composition of modules, we defined the flexibility of a node fi to be the 
number of times that a node changed modular assignment throughout the session, 
normalized by the total number of changes that were possible (i.e., by the number 
of consecutive pairs of layers in the multilayer framework). We then defined the 
flexibility of the entire network as the mean flexibility over all nodes in the network: 

See Supplementary Materials for further mathematical details and methodolog- 
ical descriptions. 
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Fig. 1. Structure of the Investigation (A) To characterize the network struc 
ture of low-frequency functional connectivity [51] at each temporal scale, we partitioned the 
raw fMRI data (left, top) from each subject's brain into signals originating from TV = 112 
cortical structures, which constitute the network's 'nodes' (right, top). The functional connec- 
tivity, constituting the network 'edges', between two cortical structures is given by a Pearson 
correlation between the mean regional activity signals (right, middle). We then statistically 
corrected the resulting TV X TV correlation matrix using a false discovery rate correction [52] to 
construct a subject-specific weighted functional brain network (left, middle). (B) Schematic of 
the investigation that was performed over the temporal scales of days, hours, and minutes. The 
complete experiment, which defines the largest scale, took place over the course of three days. 
At the intermediate scale, we conducted further investigations of the experimental sessions that 
occurred on each of those three days. Finally, to examine higher-frequency temporal structure, 
we cut each experimental session into 25 non-overlapping windows, each of which was a few 
minutes in duration. 
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Fig. 2. Multiscale Modular Architecture (A) Results for the modular decomposition of functional connectivity across temporal scales. In each panel, the 
network plots on the left show the extracted modules; different colors indicate different modules. Panels (A) and (B) correspond to the entire experiment and individual 
sessions, respectively. Boxplots show the modularity index Q (left) and the number of modules (right) in the brain network compared to randomized networks. See Methods 
for a formal definition of Q. Panel ( C) shows Q and the number of modules for the cortical (blue) compared to randomized networks (red) over the 75 time windows. Error 
bars indicate standard deviation in the mean over subjects. 
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Fig. 3. Temporal Dynamics of Modular Architecture (A) Schematic 
of a toy network with four nodes and four edges in a single time window. (B) Multilayer 
network framework in which the networks from four time windows are linked by connecting 
nodes in one time window to themselves in the adjacent time windows (colored curves). (C) 
Statistical framework composed of a connectional null model (top), a nodal null model (middle), 
and a temporal null model (bottom) in which intra-network links, inter-network links, and 
time windows respectively (shown in red) in the real network are randomized in the permuted 
network. (D) Boxplots showing differences in modular architecture between the real and 
permuted networks for the connectional (top), nodal (middle), and temporal (bottom) null 
models. We measured the structure of the network using the modularity index Q, the number 
of modules, the module size, and stationarity, which is defined as the mean similarity in the nodal 
composition of modules over consecutive time steps. Below each plot, we indicate by asterisks 
the significance of one-sample t-tests that assess whether the differences that we observed were 
significantly different from zero (gray lines): a single asterisk indicates p < .05, two asterisks 
indicate p < 1 X 10~ 6 , and three asterisks indicate p < 1 X 10 -20 . 
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Fig. 4. Flexibility and Learning (A) Boxplots showing that the increase in 
flexibility from experimental Session 1 to Session 2 was significantly greater than zero (a one- 
sample t-test gives the result t PS 6.00 with p PS 2 X 10 — 8 ), and that the magnitude 
of the decrease in flexibility from Session 2 to Session 3 was significantly greater than zero 
(t PS 7.46, p PS 2 X 10 — 1 ). f.B) Significant predictive correlations between flexibility in 
Session 1 and learning in Session 2 (black curve; p PS 0.001) and between flexibility in Session 
2 and learning in Session 3 (red curve; p PS 0.009). Note that relationships between learning 
and network flexibility in the same experimental sessions (1 and 2) were not significant; we 
obtained p > 0.13 using permutation tests. ( C) Brain regions whose flexibility in Session 1 
predicted learning in Session 2 (jp < 0.05; uncorrected for multiple comparisons). Regions that 
also passed false-positive correction were the left anterior fusiform cortex and the right inferior 
frontal gyrus, thalamus, and nucleus accumbens. (D) Brain regions whose flexibility in Session 
2 predicted learning in Session 3 (p < 0.05; uncorrected for multiple comparisons). Regions 
that also passed false-positive correction for multiple comparisons were the left intracalcarine 
cortex, paracingulate gyrus, precuneus, and lingual gyrus and the right superior frontal gyrus 
and precuneus cortex. In panels ( C-D), color indicates the Spearman correlation coefficient r 
between flexibility and learning. 
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Full Description of Methods 



Sample 

Twenty-five right-handed participants (16 female, 9 male) volunteered with informed consent in accor- 
dance with the Institutional Review Board/Human Subjects Committee, University of California, Santa 
Barbara. Handedness was determined by the Edinburgh Handedness Inventory. The mean age of the 
participants was 24.25 years (range 18.5-30 years). Of these, 2 participants were removed because their 
task accuracy was less than 60% correct, 1 was removed because of a cyst in presupplementary motor 
area (preSMA), and 4 were removed for shortened scan sessions. This left 18 participants in total. All 
participants had less than 4 years of experience with any one musical instrument, had normal vision, and 
had no history of neurological disease or psychiatric disorders. Participants were paid for their partici- 
pation. All participants completed 3 training sessions in a 5-day period, and each session was performed 
inside the Magnetic Resonance Imaging (MRI) scanner. 

Experimental Setup and Procedure 

Participants were placed in a supine position in the MRI scanner. Padding was placed under the knees in 
order to maximize comfort and provide an angled surface to position the stimulus response box. Padding 
was placed under the left forearm to minimize muscle strain when participants typed sequences. Finally, 
in order to minimize head motion, padded wedges were inserted between the participant and head coil of 
the MRI scanner. For all sessions, participants performed a cued sequence production (CSP) task (see 
Figure SI), responding to visually cued sequences by generating responses using their non-dominant (left) 
hand on a custom fiber-optic response box. For some participants, a small board was placed between the 
response box and the lap in order to help balance the box effectively. Responses were made using the 4 
fingers of the left hand (the thumb was excluded). Visual cues were presented as a series of musical notes 
on a 4-line music staff. The notes were reported in a manner that mapped the top line of the staff to the 
leftmost key depressed with the pinkie finger and so on, so that notes found on the bottom line mapped 
onto the rightmost key with the index finger (Figure SIB). Each 12-element note sequence contained 3 
notes per line, which were randomly ordered without repetition and free of regularities such as trills (e.g., 
121) and runs (e.g., 123). The number and order of sequence trials was identical for all participants. 
A trial began with the presentation of a fixation signal, which was displayed for 2 sec. The complete 
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12-element sequence was presented immediately following the removal of the fixation, and participants 
were then instructed to respond as soon as possible. They were given a period of 8 sec to type each 
sequence correctly. Participants trained on a set of 16 unique sequences, and there were three different 
levels of training exposure. Over the course of the three training sessions, three sequences — known as 
skilled sequences — were presented frequently, with 189 trials for each sequence. A second set of three 
sequences, termed familiar sequences, were presented for 30 trials each throughout training. A third set 
composed of 10 different sequences, known as novice sequences, were also presented; each novice sequence 
was presented 4-8 times during training. 

Skilled and familiar sequences were practiced in blocks of 10 trials, so that 9 out of 10 trials were 
composed of the same sequence and 1 of the trials contained a novice sequence. If a sequence was 
reported correctly, then the notes were immediately removed from the screen and replaced with the 
fixation signal, which remained on the screen until the trial duration (8 sec) was reached. If there were 
any incorrect movements, then the sequence was immediately replaced with the verbal cue INCORRECT 
and participants subsequently waited for the start of the next trial. Trials were separated with an inter- 
trial interval (ITI) lasting between sec and 20 sec, not including any time remaining from the previous 
trial. Following the completion of each block, feedback (lasting 12 sec and serving as a rest) was presented 
that detailed the number of correct trials and the mean time that was taken to complete a sequence. 
Training epochs contained 40 trials (i.e., 4 blocks) and lasted a total of 345 scan repetition times (TRs), 
which took a total of 690 sec. There were 6 scan epochs per training session (2070 scan TRs). In total, 
each skilled sequence was presented 189 times over the course of training (18 scan epochs; 6210 TRs). 

In order to familiarize participants with the task, they were given a short series of warm-up trials 
the day before the initial training session inside the scanner. Practice was also given in the scanner 
during the acquisition of the structural scans and just prior to the start of the first training-session 
epoch. Stimulus presentation was controlled with MATLAB® version 7.1 (Mathworks, Natick, MA) in 
conjunction with Cogent 2000 (Functional Imaging Laboratory, 2000). Key-press responses and response 
times were collected using a fiber-optic custom button box transducer that was connected to a digital 
response card (DAQCard-6024e; National Instruments, Austin, TX). We assessed learning using the slope 
of the movement time (MT), which is the difference between the time of the first button press and the 
time of the last button press in a single sequence (see Figure SIB) [1]. The negative slope of the movement 
curve over trials indicates that learning is occurring [lj. 
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Acquisition and Preprocessing of fMRI Data 

Functional MRI (fMRI) recordings were conducted using a 3.0 T Siemens Trio with a 12-channel phased- 
array head coil. For each functional run, a single-shot echo planar imaging that is sensitive to blood 
oxygen level dependant (BOLD) contrast was used to acquire 33 slices (3 mm thickness) per repetition 
time (TR), with a TR of 2000 ms, an echo time (TE) of 30 ms, a flip angle of 90 degrees, and a field of 
view (FOV) of 192 mm. The spatial resolution of the data was defined by a 64 x 64 acquisition matrix. 
Before the collection of the first functional epoch, a high-resolution Tl-weighted sagittal sequence image 
of the entire brain was acquired (TR = 15.0 ms, TE = 4.2 ms, flip angle = 9 degrees, 3D acquisition, 
FOV = 256 mm; slice thickness = 0.89 mm, and spatial acquisition matrix dimensions = 256 x 256). 

All image preprocessing was performed using the FMRIB (Oxford Centre for Functional Magnetic 
Resonance Imaging of the Brain) Software Library (FSL) [2] . Motion correction was performed using the 
program MCFLIRT (Motion Correction using FMRIB 's Linear Image Registration Tool). Images were 
high-pass filtered with a 50 sec cutoff period. Spatial smoothing was performed using a kernel where 
the full width at half maximum was 8 mm. No temporal smoothing was performed. The signals were 
normalized globally to account for transient fluctuations in signal intensity. 

Partitioning the Brain into Regions of Interest 

Brain function is characterized by a spatial specificity: different portions of the cortex emit inherently 
different activity patterns that depend on the experimental task at hand. In order to measure the 
functional connectivity between these different portions, it is common to apply an atlas of the entire 
brain to raw fMRI data in order to combine information from all 3 mm cubic voxels found in a given 
functionally or anatomically defined region (for recent reviews, see [3-5]). Several atlases are currently 
available, and each provides slightly different parcellations of the cortex into discrete volumes of interest. 
Several recent studies have highlighted the difficulty of comparing results from network analyses derived 
from different atlases [6-8]. In the present work, we have therefore used a single atlas that provides the 
largest number of uniquely identifiable regions — this is the Harvard-Oxford (HO) atlas, which is available 
through the FSL toolbox [2,9]. The HO atlas provides 112 functionally and anatomically defined cortical 
and subcortical regions; for a list of the brain regions, see Supplementary Table 1. Therefore, for each 
individual fMRI data set, we estimated regional mean BOLD time series by averaging voxel time series 
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in each of the 112 regions. Each regional mean time series was composed of 2070 time points for each of 
the 3 experimental sessions (for a total of 6210 time points for the complete experiment). 

Wavelet Decomposition 

Brain function is also characterized by a frequency specificity; different cognitive and physiological func- 
tions are associated with different frequency bands, which can be investigated using wavelets. Wavelet 
decompositions of fMRI time series have been applied extensively in both resting-state and task-based con- 
ditions [10, 11]. In both cases, they provide increased sensitivity for the detection of small signal changes 
in non-stationary time series with noisy backgrounds [12] . In particular, the maximum-overlap discrete 
wavelet transform (MODWT) has been extensively used in connectivity investigations of fMRI [13-18]. 
Accordingly, we used MODWT to decompose each regional time series into wavelet scales corresponding 
to specific frequency bands [19]. We were interested in quantifying high-frequency components of the 
fMRI signal, correlations between which might be indicative of cooperative temporal dynamics of brain 
activity during a task. Because our sampling frequency was 2 sec (1 TR = 2 sec), wavelet scale one 
provided information on the frequency band 0.125-0.25 Hz and wavelet scale two provided information 
on the frequency band 0.06-0.125 Hz. Previous work has indicated that functional associations between 
low- frequency components of the fMRI signal (0-0.15 Hz) can be attributed to task-related functional 
connectivity, whereas associations between high-frequency components (0.2-0.4 Hz) cannot [20]. This 
frequency specificity of task-relevant functional connectivity is likely to be due at least in part to the 
hemodynamic response function, which might act as a noninvertible bandpass filter on underlying neural 
activity [20]. In the present study, we therefore restricted our attention to wavelet scale two in order to 
assess dynamic changes in task-related functional brain architecture over short time scales while retaining 
sensitivity to task-perturbed endogenous activity [21], which is most salient at about 0.1 Hz [22-24]. 

Connectivity Over Multiple Temporal Scales 

Multiscale Connectivity Estimation We measured functional connectivity over three temporal 
scales: the large scale of the complete experiment (which lasted 3 hours and 27 minutes), the session 
time scale of each fMRI recording session (3 sessions of 69 minutes each; each session corresponded to 
2070 time points), and the shorter time scales of intra-session time windows (where each time window 
was approximately 3.5 min long and lasted 80 time points). 
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In the investigation of large-scale connectivity, we concatenated regional mean time series over all 3 
sessions, as has been done previously [25]. We then constructed for each subject a functional association 
matrix based on correlations between regional mean time series. At the mesoscopic scale, we extracted 
regional mean time series from each experimental session separately to compute session-specific matrices. 
At the small scale, we constructed intra-session time windows with a length of T = 80 time points, 
giving a total of 25 time windows in each session (see the Results section of this supplementary document 
for a detailed investigation across a range of T values). We constructed separate functional association 
matrices for each subject in each time window (25) for each session (3) for a total of 75 matrices per 
subject. We chose the length of the time window to be long enough to allow adequate estimation of 
correlations over the frequencies that are present in the wavelet band of interest (0.06-0.12 Hz), yet short 
enough to allow a fine-grained measurement of temporal evolution over the full experiment. 

Construction of Brain Networks To construct a functional network, we must first define a measure 
of functional association between regions. Measures of functional association range from simple linear 
correlation to nonlinear measures such as mutual information. In the majority of network investigations in 
fMRI studies to date, the measure of choice has been the Pearson correlation [13,15,18,26,27], perhaps due 
to its simplicity and ease of interpretation. Therefore, in order to estimate static functional association, 
we calculated the Pearson correlation between the regional mean time series of all possible pairs of regions 
i and j. This yields an N x N correlation matrix with elements r^-, where N = 112 is the number of 
brain regions of interest in the full brain atlas (see earlier section on "Partitioning the Brain into Regions 
of Interest" for further details). 

However, as pointed out in other network studies of fMRI data [13], not all elements r it j of the full 
correlation matrix necessarily indicate significant functional relationships. Therefore, in addition to the 
correlation matrix element nj, we computed the p- value matrix element pij, which give the probabilities 
of obtaining a correlation as large as the observed value r^j by random chance when the true correlation 
is zero. We estimated p-values using approximations based on the ^-statistic using the MATLAB® 
function corrcoef [28]. In the spirit of Ref. [29] and following Ref. [13], we then tested the p- values p it j 
for significance using a False Discovery Rate (FDR) of p < 0.05 to correct for multiple comparisons [30,31] . 
We retained matrix elements j whose p-values ptj passed the statistical FDR threshold. Elements of 
rij whose p- values pij did not pass the FDR threshold were set to zero in order to create new correlation 
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matrix elements r' „•. 

We applied the statistical threshold to all r^j independent of the sign of the correlation. Therefore, 
the resulting r\ ■ could contain both positive and negative elements if there existed both positive and 
negative elements of r^j whose p-values pij passed the FDR threshold. Because this was a statistical 
threshold, the network density of (defined as the fraction of non-zero matrix elements) was determined 
statistically rather than being set a priori. Network density varied over temporal resolutions; the mean 
density and standard deviation for networks derived from correlation matrices at the largest time scale 
(3 hr and 27 minutes) was 0.906 (0.019%), at the intermediate time scale (69 min) was 0.846 (0.029), and 
at the short time scale (3.5 min) was 0.423 (0.110). 

We performed the procedure described above for each subject separately to create subject-specific 
corrected correlation matrices. These statistically corrected matrices gave adjacency matrices A (see the 
discussion below) whose elements were Aij = ■. 

Network Modularity To characterize the large-scale functional organization of the subject-specific 
weighted matrices A, we used tools from network science [32]. In a network framework, brain regions 
constitute the nodes of the network, and inter-regional functional connections that remain in the connec- 
tivity matrix constitute the edges of the network. One powerful concept in the study of networks is that 
of community structure, which can be studied using algorithmic methods [33,34]. Community detection 
is an attempt to decompose a system into subsystems (called 'modules' or 'communities'). Intuitively, a 
module consists of a group of nodes (in our case, brain regions) that are more connected to one another 
than they are to nodes in other modules. A popular way to investigate community structure is to optimize 
the partitioning of nodes into modules such that the quality function Q is maximized (see [33, 34] for 
recent reviews and [35] for a discussion of caveats), for which we give a formula below. 

From a mathematical perspective, the quality function Q is simple to define. One begins with a graph 
composed of N nodes and some set of connections between those nodes. The adjacency matrix A is then 
an N x TV matrix whose elements Ay detail a direct connection or 'edge' between nodes i and j, with a 
weight indicating the strength of that connection. The quality of a hard partition of A into communities 
(whereby each node is assigned to exactly one community) is then quantified using the quality function 
Q. Suppose that node i is assigned to community gi and node j is assigned to community gj. The most 
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popular form of the quality function takes the form [33, 34] 



(i) 



where S(gi,gj) = 1 if = <? 7 and it equals otherwise, and P^j is the expected weight of the edge 
connecting node i and node j under a specified null model. (The specific choice of Q in Equation 1 is 
called the network modularity or modularity index [36].) A most common null model (by far) used for 
static network community detection is given by [33,34,37] 



where ki is the strength of node i, kj is the strength of node j, and m — \ Ay. The maximization 
of the modularity index Q gives a partition of the network into modules such that the total edge weight 
inside of modules is as large as possible (relative to the null model, subject to the limitations of the 
employed computational heuristics, as optimizing Q is NP-hard [33,34,38]). 

Network modularity has been used recently for investigations of resting-state functional brain networks 
derived from fMRI [26,27] and of anatomical brain networks derived from morphometric analyzes [39]. In 
these previous studies, brain networks were constructed as undirected binary graphs, so that each edge 
had a weight of either 1 or 0. The characteristics of binary graphs derived from neuroimaging data are 
sensitive to a wide variety of cognitive, neuropsychological, and neurophysiological factors [4,5]. However, 
increased sensitivity is arguably more likely in the context of the weighted graphs that we consider, as 
they preserve the information regarding the strength of functional associations (though, as discussed 
previously, matrix elements rjj that are statistically insignificant are still set to 0) [40]. An additional 
contrast between previous studies and the present one is that (to our knowledge) investigation of network 
modularity has not yet been applied to task-based fMRI experiments, in which modules might have a 
direct relationship with goal-directed function. 

We partitioned the networks represented by the weighted connectivity matrices into n communities by 
using a Louvain greedy community detection method [41] to optimize the modularity index Q. Because 
the edge weights in the correlation networks that we constructed contain both positive and negative 
correlation coefficients, we used the signed null model proposed in Ref. [42] to account for communities of 



k^ kj 

2^r ' 



(2) 
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nodes associated with one another through both negative and positive edge weights. (Recall that we are 
presently discussing aggregated correlation networks A, so we are detecting communities in single-layer 
networks, as has been done in previous work. In order to investigate time-evolving communities, we will 
later employ a new mathematical development that makes it possible to perform community detection in 
multilayer networks [43].) We first defined w f- to be an TV x N matrix containing the positive elements 
of Aij and w~- to be an N x TV matrix containing only the negative elements of Aij. The quality function 
to be maximized is then given by 



where gi is the community to which node i is assigned, gj is the community to which node j is assigned, 



resolution parameter values to unity. 

In our investigation, we have focused on the mean properties of ensembles of partitions rather than 
on detailed properties of individual partitions. This approach is consistent with recent work illustrating 
the fact that the optimization of quality functions like Q and Q± is hampered by the complicated shape 
of the optimization landscape. In particular, one expects to find a large number of partitions with near- 
optimum values of the quality function [35], collectively forming a high-modularity plateau. Theoretical 
work estimates that the number of "good" (in the sense of high values of Q and similar quality functions) 
partitions scales as 2 n_1 , where n is the mean number of modules in a given partition [35]. In both toy 
networks and networks constructed from empirical data, many of the partitions found by maximizing 
a quality function disagree with one another on the components of even the largest module, impeding 
interpretations of particular partitions of a network [35] . Therefore, in the present work, we have focused 
on quantifying mean qualities of the partitions after extensive sampling of the high-modularity plateau. 
Importantly, the issue of extreme near-degeneracy of quality functions like Q is expected to be much less 
severe in the networks that we consider than is usually the case, because we are examining small, weighted 
networks rather than large, unweighted networks [35]. We further investigate the degenerate solutions in 
terms of their mean, standard deviation, and maximum. We find that Q± values are tightly distributed, 
with maximum values usually less than three standard deviations from the mean (see Supplementary 
Results). 




(3) 



7 + and 7 are resolution parameters, and = J2j w tj> w i = 



J2j w ij [42] . For simplicity, we set the 
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Statistical Testing To determine whether the value of Q prn or the number of modules was greater 
or less than expected in a random system, we constructed randomized networks with the same degree 
distribution as the true brain networks. As has been done previously [27,44], we began with a real brain 
network and then iteratively rewired it using the rewiring algorithm of Maslov and Sneppen [45] . The 
procedure we used for accomplishing this rewiring was to choose at random two edges — one that connects 
node A to node B and another that connects nodes C and D — and then to rewire them to connect A 
to C and B to D. This allows us to preserve the degree, or number of edges, emanating from each node 
although it does not retain a node's strength. To ensure a thorough randomization of the underlying 
connectivity structure, we performed this procedure multiple times, such that the expected number of 
times that each edge was 'rewired' was 20. This null model will be hereafter referred to as the static 
random network null model. (This is distinct from the null models that we have developed for statistical 
testing of community structure in multilayer networks, as discussed in the main manuscript and in later 
sections of this Supplement.) The motivation for this process is to compare the brain with a null model 
that resembles the configuration model [46] , which is a random graph with prescribed degree distribution. 

We constructed 100 instantiations of the static random network null model for each real network that 
we studied. We constructed representative values for diagnostics from the random networks by taking 
the mean network modularity and mean number of modules over those 100 random networks. We then 
computed the difference between the representative random values and the real values for each diagnostic, 
and we performed a one-sample i-test over subjects to determine whether that difference was significantly 
greater than or less than zero. For each case, we then reported p-values for these tests. 

Sampling of the static random network null model distribution is important in light of the known 
degeneracies of modularity (which we discuss further in the Supplementary Results section below) [35]. 
One factor that accounts for a significant amount of variation in Q± is the size (i.e., number of nodes) 
of the network, so comparisons between networks of different sizes must be performed with caution. 
Therefore, we note that all networks derived from the aforementioned null model retain both the same 
number of nodes and the same number of edges as the real networks under study. This constrains 
important factors in the estimation of Q± . 

Visualization of Networks We visualized networks using the software package MATLAB® (2007a, 
The Math Works Inc., Natick, MA). Following Ref. [47], we used the Fruchtermann-Reingold algorithm [48] 
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to determine node placement for a given network with respect to the extracted communities and then 
used the Kamada-Kawai algorithm [49] to place the nodes within each community. 

Multilayer Network Modularity: Temporal Dynamics of Intra-Session Con- 
nectivity 

In order to investigate the temporal evolution of modular architecture in human functional connectivity, 
we used a mutilayer network framework in which each layer consists of a network derived from a single 
time window. Networks in consecutive layers therefore correspond to consecutive time windows. We 
linked networks in consecutive time windows by connecting each node in one window to itself in the 
previous and in the next windows (as shown in Figure 3A-B in the main text) [43]. We constructed a 
multilayer network for each individual and in each of the three experimental sessions. We then performed 
community detection by optimizing a multilayer modularity (see the discussion below) [43] using the 
Louvain greedy algorithm (suitably adapted for this more general structure) on each multilayer network 
in order to assess the modular architecture in the temporal domain. 

Our examination of static network architecture, we used the wavelet correlation to assess functional 
connectivity. Unfortunately, more sensitive measures of temporal association such as the spectral co- 
herence are not appropriate over the long time scales assessed in the static investigation due to the 
nonstationarity of the fMRI time series [10-12], and it is exactly for this reason that we have used the 
wavelet correlation for the investigation of aggregated (static) networks. However, over short temporal 
scales such as those being used to construct the multilayer networks, fMRI signals in the context of the 
motor learning task that we study can be assumed to be stationary [50] , so spectral measures such as the 
coherence are potential candidates for the measurement of functional association. 

In the examination of the dynamic network architecture of brain function using multilayer community 
detection, our goal was to measure temporal adaptivity of modular function over short temporal scales. 
In order to estimate that temporal adaptivity with enhanced precision, we used the magnitude squared 
spectral coherence (as estimated using the minimum- variance distortionless response method [51]) as a 
measure of nonlinear functional association between any two time series. In using the coherence, which 
has been demonstrated to be useful in the context of fMRI neuroimaging data [20], we were able to 
measure frequency-specific linear relationships between time series. 

As in the static network analysis described earlier, we tested the elements of each N x N coherence 
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matrix (which constitutes a single layer) for significance using an FDR correction for multiple comparisons. 
We used the original weighted (coherence values) of network links corresponding to the elements that 
passed this statistical test, while those corresponding to elements that did not pass the test were set to 
zero. In applying a community detection technique to the resulting coherence matrices, it is important to 
note that the coherence is bounded between and 1. We can therefore use a multilayer quality function 
with an unsigned null model rather than the signed null model used in the static case described earlier. 
The multilayer modularity Q m i is given by [43] 



Qmi = y ^2 I \ - H-^^-j s ir + 5 lJ C j i r \ 8(gu,g jr 

ijlr 



), (4) 



where the adjacency matrix of layer I (i.e., time window number I) has components Am, 7; is the resolution 
parameter of layer I, gu gives the community assignment of node i in layer I, gj r gives the community 
assignment of node j in layer r, Cji r is the connection strength between node j in layer r and node j in 
layer I (see the discussion below), ku is the strength of node i in layer I, 2ji = Ylj r K jr, K ji — kji+Cji, and 
cji — J2 r Cjir- For simplicity, as in the static network case, we set the resolution parameter 7; to unity 
and we have set all non-zero Cji r to a constant C, which we will term the 'inter-layer coupling'. In the 
main manuscript, we report results for C = 1. In the Supplementary Results section of this document, 
we investigate the dependence of our results on alternative choices for the value of C. 

Diagnostics We used several diagnostics to characterize dynamic modular structure. These include 
the multilayer network modularity Q m u the number of modules n, the module size s, and the stationarity 
of modules (. We defined the size of a module s to be the mean number of nodes per module over all time 
windows over which the community exists. We used the definition of module stationarity from Ref. [52]. 
We started by calculating the autocorrelation function U(t) of two states of the same community G(t) at 
t = 1 time steps apart using the formula 

= \G(t )nG(t + t)\ 
U[t) ~ \G(t )UG(t +t)\ ' [b) 

where t is the time at which the community is born, \G(t Q ) n G(t n +t)\ is the number of nodes that are 
members of both G(to) and G(to + t), and \G(t ) U G(to +t)\ is the total number of nodes in the union 
of G(to) and G(to +t) [52]. We defined t' to be the final time step before the community is extinguished. 
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The stationarity of a community is then given by 

Et" t p u(t,t + i) 

^ t'-to-l ' 

which is the mean autocorrelation over consecutive time steps [52]. 

Statistical Framework The study of the "modular architecture" of a system is of little value if the 
system is not modular. It is therefore imperative to statistically quantify the presence or absence of 
modular architecture to justify the use of community detection in a given application. Appropriate 
random null models have been developed and applied to the static network framework [27, 44] , but no 
such null models yet exist for the multilayer framework. We therefore developed several null models in 
order to statistically test the temporal evolution of modular structure. We constructed three independent 
null models to test for (1) network structure dependent on the topological architecture of connectivity, 
(2) network structure dependent on nodal identity, and (3) network structure dependent on the temporal 
organization of layers in the multilayer framework. 

In the connectional null model (1), we scrambled links between nodes in any given time window (the 
entire experiment, 3.45 hr; the individual scanning session, 69 min; or intra-session time windows, 3.45 
min) while maintaining the total number of connections emanating from each node in the system. To 
be more precise, for each layer of the multilayer network, we sampled the static random network null 
model (see the discussion above in the context of static connectivity architecture) for that particular 
layer. That is, we reshuffled the connections within each layer separately while maintaining the original 
degree distribution. We then linked these connectivity-randomized layers together by coupling a node in 
one layer to itself in contiguous layers to create the connectional null model multilayer network, just as we 
connected the real layers to create the real multilayer network. In the present time-dependent context, 
we performed this procedure on each time window in the multilayer network, after which we applied 
the multilayer community detection algorithm to determine the network modularity of the randomized 
system. 

In constructing a nodal null model (2), we focused on the links that connected a single node in one 
layer of the multilayer framework to itself in the next and previous layers. In the null model, the links 
between layers connect a node in one layer to randomly-chosen nodes in contiguous layers instead of 
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(6) 



connecting the node to itself in those layers. Specifically, in each time window Tj (except for the final 
one), we randomly connected the nodes in the corresponding layer to other nodes in the next time window 
(Tj + i) such that no node in Tj was connected to more than one node in r i+ i. We then connected nodes 
in Tj+i to randomly-chosen nodes in Tj+2, and so on until links between all time windows had been fully 
randomized. 

We also considered randomization of the order in which time windows were placed in the multilayer 
network to construct a temporal null model (3). In the real multilayer construction, we (of course) 
always placed the network from time window just before the network from time window Tj+i. In 
the temporal null model, we randomly permuted the temporal location of the individual layer in the 
multilayer framework such that the probability of any time window Tj following any other time window 
r i U i) was uniform. 

Statistical Testing In both the real network and the networks derived from null models (l)-(3), it is 
important to adequately sample the distributions of partitions meant to optimize the modularity index 
Qmi . This step in our investigation was particularly important in light of the extreme degeneracy of the 
network modularity function Q m i [35] (see the Supplementary Results section on the Degeneracies of Q 
for a quantitative characterization of such degeneracies). 

Because the multilayer community detection algorithm can find different maxima each time it is run, 
we computed the community structure of each individual real multilayer network a total of 100 times. 
We then averaged the values of all diagnostics (modularity index Q m i, number of modules, module size, 
and stationarity) over those 100 partitions to create a representative real value. To perform our sampling 
for the null models, we considered 100 multilayer network instantiations for each of the three different 
null models. We also performed community detection on these null models using our multilayer network 
adaptation of the Louvain modularity-optimization algorithm [41] to create a distribution of values for 
each diagnostic. We then used the mean value of each diagnostic in our subsequent investigation as the 
representative value of the null model. 

We used one-sample t-tests to test statistically whether the differences between representative values 
from the real networks and null model networks over the subject population was significantly different 
from zero. The results, which we reported in the main manuscript, indicated that in contrast to what we 
observed using each of the three null models, the human brain displayed a heightened modular structure. 
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That is, it is composed of more modules, which have smaller sizes. Considering the three null models in 
order, this suggests that cortical connectivity has a precise topological organization, that cortical regions 
consistently maintain individual connectivity signatures necessary for cohesive community organization, 
and that functional communities evolve cohesively in time (see Figure 2 in the main manuscript). Im- 
portantly, the stationarity of modular organization £ was also higher in the human brain than in the 
connectional or nodal null models, indicating a cohesive temporal evolution of functional communities. 

Temporal Dynamics of Brain Architecture and Learning 

In the present study, we have attempted to determine whether changes in the dynamic modular architec- 
ture of functional connectivity is shaped by learning. We assessed the learning in each session using the 
slope of the movement times (MT) of that session. Movement time is defined as the difference between 
the time of the first button press and the time of the last button press in a single sequence (see Figure 
SIB). During successful learning, movement time is known to fall logarithmically with time [1]. However, 
two subjects from session 1 and one subject from session 2 showed an increasing movement time as the 
session progressed. We therefore excluded these three data points in subsequent comparisons due to the 
decreased likelihood that successful learning was taking place. This process of screening participants 
based on movement time slope is consistent with previous work suggesting that fMRI activation patterns 
during successful performance might be inherently different when performance is unsuccessful [53]. 

In principle, modular architecture might vary with learning by displaying changes in global diagnostics 
such as the number of modules or the modularity index Q or by displaying more specific changes in the 
composition of modules. To measure changes in the composition of modules, we defined the flexibility 
of a node /; to be the number of times that node changed modular assignment throughout the session, 
normalized by the total number of changes that were possible (i.e., by the number of consecutive pairs 
of layers in the multilayer framework). We then defined the flexibility of the entire network as the mean 
flexibility over all nodes in the network: F = Y^iLi Si- 

Statistics and Software 

We implemented all computational and simple statistical operations using the software packages MATLAB® 
(2007a, The MathWorks Inc., Natick, MA) and Statistica® (version 9, StatSoft Inc.). We performed the 
network calculations using a combination of in-house software (including multilayer community detection 
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code [43]) and the Brain Connectivity Toolbox [40]. 
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Supplementary Results 



Degeneracies of Q 

As discussed earlier in the Methods section, we focused in this investigation on the mean properties of 
ensembles of partitions rather than on detailed properties of individual partitions. Our approach was 
motivated by recent work indicating that the optimization of modularity and similar quality functions 
is hampered by the complicated shape of the optimization landscape, which includes a large number of 
partitions with near-optimum values that collectively form a high modularity plateau [35]. To quantify 
and address this degeneracy of Q± and Q m i, we now provide supplementary results on the mean, standard 
deviation, and maximum values of Q± and Q m i over the 100 samples of the plateau computed for all real 
networks in both the static and dynamic frameworks. 

The mean number of modules in a given partition in the static framework was n « 3.08 for the entire 
experiment, n « 3.07 for individual experimental sessions, and n « 3.55 for the small intra-session time 
windows. The mean number of modules in a given partition in the multilayer framework was n « 6.00. 
We have therefore chosen to sample the quality functions Q± and Q m i a total of 100 times (which is more 
than 2™ _1 in each case, and therefore adequately samples the degenerate near-optimum values of Q± and 
Qmi [35]). In order to characterize the distribution of solutions found in these 100 samplings, we have 
computed the mean, standard deviation, and maximum of Q± (static cases) and Q m i (dynamic cases); 
see Figure S2. We found that the values of Q± and Q m i are tightly distributed, and that the maximum 
values of Q± or Q m i are between and 3 standard deviations higher than the mean. Although we 
remain cautious because we have not explored all possible computational heuristics, we are nevertheless 
encouraged by these results that the mean values of Q± and Q m i that we have reported are representative 
of the true maximization of the two quality functions. 

Reproducibility We calculated the intra-class correlation coefficient (ICC), to determine whether 
values of Q± and Q m i derived from a single individual over the 100 samples were more similar than 
values of Q± or Q m i derived from different individuals. The ICC is a measure of the total variance for 
which between-subject variation accounts [54, 55] , and it is defined as 



ICC = 



.2 

bs 



(7) 
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where at, s is the between-subject variance and a ws is the pooled within-subject variance ('pooled' indicates 
that variance was estimated for each subject and then averaged over subjects). The ICC is normalized 
to have a maximum value of 1; values above 0.5 indicate that there is more variability between Q± and 
Q m i values from different subjects than between Q± and Q m i values from the same subject. In the static 
framework, the ICC was 0.9884 at the large scale (the entire experiment), an average of 0.9863 at the 
intermediate scale (three experimental sessions), and an average of 0.9847 at the small scale (individual 
time windows). In the multilayer framework, we calculated that ICC « 0.9983. These results collectively 
indicate that the Q± and Q m i values that we reported in this work were significantly reproducible over 
the 100 samples of the respective quality function landscape. That is, the Q± or Q m i values drawn from 
the 100 samples of a single subject's network modularity landscape were more similar than Q± or Q m i 
values drawn from different subjects. 

Effect of the Inter-Layer Coupling Parameter 

The multilayer network framework requires one to define a coupling parameter C that indicates the 
strength of the connections from a node in one time window to itself in the two neighboring time windows 
[43] . In order to be sensitive to both temporal dynamics and intra- layer network architecture, the coupling 
parameter should be on the same scale of values as the edge weights. For example, if edge weights are 
coherence values lying between and 1, then the coupling parameter also ought to lie between and 1. 
In the results that we presented in the main manuscript, we set the coupling parameter to be C = 1, 
which is the highest value consistent with the intra-layer edge weights given by the normalized coherence. 
However, if we were to alter the coupling value, one might expect the number of communities to be 
altered in kind. As the strength of the coupling is increased, one might expect fewer communities to be 
uncovered due to the increased temporal dependence between layers [43]. Similarly, as the inter-layer 
coupling is weakened, one might expect more communities to be detected. 

To probe the effect of the inter-layer coupling strength, we thus varied C from well below to well above 
the maximum intra-layer edge weight (0.2 < C < 2). In Figure S3 (cortical network results are shown 
in blue), we illustrate the effects of sweeping over this coupling parameter on our four diagnostics. The 
modularity index Q m ; increases with increasing inter-layer coupling, whereas the other three diagnostics — 
number of modules, module size, and stationarity — increase initially and then plateau approximately at 
about C = 1 and above. The change in behavior near C = 1 can be rationalized as follows: For C < 1, 
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intra-layer edge weights dominate the modularity optimization, whereas inter-layer edge weights dominate 
for C > 1. The proposed choice of C = 1 therefore balances the impact of known coherence in brain 
activity (as given by the intra-layer edge weights) on measured architectural adaptations and is therefore 
a natural choice with which to investigate biologically meaningful organization. 

We also computed 100 temporal, nodal, and connectional null model networks for each of the additional 
coupling parameter values (see Figure S3; null model network results shown in green, orange, and red). 
The results indicate that the relationship between diagnostics in the cortical networks and null model 
networks is dependent on the diagnostic. For example, modularity values of null model networks are 
consistently lower than modularity values of cortical networks. However, stationarity in the null model 
networks is lower than that in cortical networks for small values of G but higher than that of cortical 
networks for high values of C. This nontrivial behavior suggests an added sensitivity of the proposed 
null model networks to the multilayer network construction, which might be useful in other experimental 
contexts and therefore warrants further investigation. 

Effect of the Time Window Length 

In the construction of networks at the smallest time scale, it is necessary to choose a length of the time 
window T. In choosing this time window length, two considerations are important: (1) the time window 
must be short enough to adequately measure temporal evolution of network structure, and (2) the time 
window must be long enough to adequately estimate the functional association between two time series 
using (for example) the correlation or coherence [56]. In the main text, we reported results for time 
windows of 80 data points in length. This gives 25 time windows in each experimental session, for a 
total of 75 time windows over the 3 sessions. In addition to this extensive coverage of the underlying 
temporal dynamics, the choice of a time window of 80 data points in length also ensures that 20 data 
points can be used for the estimation of the functional association between time series in the frequency 
band of interest — i.e., at wavelet scale two (0.06-0.12 Hz). If one were to increase the time window 
length, one would expect a decreased ability to measure temporal variations due to the presence of fewer 
time windows per session. If one were to decrease the time window length, one would expect increased 
variance in the estimation of the functional association between time series due to the use of fewer data 
points in the estimation of either the coherence or the correlation [16]. 

To probe the effect of the time window length, we varied T from T = 80 to T = 110 (see Figure 
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S4; cortical network results are shown in blue). We find that the stationarity of the modules increases 
with increasing time window length. As T is increased, the functional association between any two nodes 
is averaged over a longer time series, so small adaptations over shorter time scales can no longer be 
measured. This smoothing is likely the cause of the increased stationarity that we find at high values of 
T. It suggests that functional association measured over long time windows is less dependent on the time 
window being used than functional association measured over short time windows. This finding supports 
our choice of short time windows in order to measure dynamic adaptations in network architecture. 

We also computed 100 temporal, nodal, and connectional null model networks for each of the additional 
time window lengths (see Figure S4; null model network results are shown in green, orange, and red). 
The results indicate that the relationships between diagnostics in the cortical networks and null model 
networks are largely conserved across time window lengths. 

Learning and Flexibility 

In the main text, we reported a significant correlation between the flexibility of dynamic modular archi- 
tecture in a given experimental session, as measured by the (normalized) number of times a node changes 
module allegiance, and learning in the subsequent experimental sessions, as measured by the slope of 
the movement time (see Methods). We found that the mean value of flexibility was approximately 0.30, 
that it fluctuated over the three experimental sessions, and that the values were highest in the second 
experimental session (see Table 2 in this Supplement). We followed this large-scale calculation with an 
investigation into the relationship between nodal flexibility (in particular brain regions) and learning. We 
found, as shown Figure 4 of the main manuscript, that the flexibility of a large number of brain regions 
could be used to predict learning in the following session. Here we also note that these regions were not 
those with highest flexibility or lowest flexibility in the brain. In fact, the flexibility of those regions that 
predicted learning was not significantly different from the flexibility of those regions that did not predict 
learning: t « 0.01 p « 0.98 (Session 1) and t « 0.87, p « 0.38 (Session 2). 

In addition to those results reported in the main manuscript, we tested whether the flexibility of 
the cortical networks was significantly different from the flexibility expected in the (connectional, nodel, 
and temporal) random network null models. As we show in Table 3 in this Supplement, the flexibility 
of the connectional and nodal null model networks was significantly higher than that of the cortical 
networks, and we found no discernible differences between the cortical networks and the temporal null 
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model networks. We found the greatest degree of flexibility in the nodal null model, in which individual 
nodes in any given time window were coupled to randomly selected nodes in the following time window. 
It is thus plausible that the subsequent disruption of nodal identity caused nodes to change computed 
module allegiances in this null model. 

Robustness to Alternative Definitions It is important to assess the robustness of our findings 
to different definitions of flexibility. We therefore defined an alternative flexibility measure f- to be 
the number of communities (modules) to which a node belongs at some point in a given experimental 
session. The mean alternative flexibility F' is then given by averaging f[ over all nodes in the network: 
F' = y^!j—i fi- Using this alternative definition of flexibility, we again tested for differences between the 
cortical network and the three random network null models. As shown in Tables 2 and 3, the F 1 values 
of cortical networks were also significantly different from those in the null model networks. Interestingly, 
for this alternative definition, the temporal network null model exhibits significantly lower flexibility 
than the cortical networks, suggesting that this measure of flexibility might be sensitive to biologically 
relevant temporal evolution of modular architecture. Finally, we tested whether this alternative definition 
of flexibility also displayed a relationship to learning. Flexibility and learning were not significantly 
correlated in Session l(r« 0.02, p w 0.90) or in Session 2 (r « 0.18, p « 0.48), but flexibility in Session 1 
was predictive of learning in Session 2 (r « 0.64, p « 0.002), and flexibility in Session 2 was predictive of 
learning in Session 3 (r « 0.51, p « 0.019). These results for the alternative flexibility F' are consistent 
with those of the original definition F, suggesting that our findings are robust. 
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Supplementary Discussion 



Resolution Limit of Modularity 

When detecting communities by optimizing modularity and similar quality functions, it is important to 
note that modularity suffers from a resolution limit [33-35,57]. As a result, the maximum-modularity 
partition can be biased towards a particular module size and can have difficulty resolving modules smaller 
than that size. Consequently, small modules of potential interest have the potential to be hidden within 
larger groups of nodes that have been detected. Modularity's resolution limit is particularly prevalent in 
sparse networks, binary networks, and large networks, and its effects tend to be much less significant in 
networks of the type (dense, weighted, and small) that we have studied [35]. 

Measuring Differences in Brain States 

In the present work, we have characterized differences in brain states during learning by examining the 
global network architecture and measuring differences between that architecture over three experimental 
sessions. An alternative line of investigation would be to seek network motifs (i.e., small patterns of nodes 
and edges) that have the potential to distinguish between brain states. This could be done using statistical 
methods [58], machine-learning techniques [59], or a combination of the two [60]. Our approach, however, 
has the advantage of assessing alterations in large-scale achitectural properties rather than differences 
in small parts of that architecture. Additionally, the approach that we have chosen provides a direct 
characterization of the underlying functional connectivity architecture irrespective of differences between 
brain states. Using this approach, we have therefore been able to demonstrate, for example, that there 
is significant non-random modular organization across multiple temporal scales. 

A Note on Computational Time 

The investigations that we reported in the present work involved about f 0, 000 CPU-days, and our 
study was therefore made possible by the use of two computing clusters available at the Institute for 
Collaborative Biotechnologies at UC Santa Barbara. Cluster f was composed of 42 Dell SCf425s (dual 
single-core Xeon 2.8GHz, 4GB memory), 5 Dell PEf950s (dual quad-core Xeon E5335 2.0GHz, 8GB 
memory), f Dell 2850 (RAID storage includes 500GB for the home directory), and MATLAB® MDCE 
with 128 worker licenses (cluster currently has 124 compute cores), Gigabit Ethernet, Software RAID 
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backup node (converted compute node) with 673GB software RAID backup. Cluster 2 was composed of 
20 HP Proliant DL160 G6s (dual quad-core E5540 "Nehalem" 2.53GHz, 24GB memory), 1 HP DL180 
G6 (RAID storage includes 2.1TB for the home directory), MATLAB® MDCE with 160 worker licenses 
(cluster currently has 160 compute cores), Gigabit Ethernet, and a storage node with 4.6TB of RAID 
storage (for backup). 

We performed maximization of the quality functions (Q pm , Qmi) a total of 100 times for every 
connectivity matrix under study. In the static connectivity investigation, we constructed connectivity 
matrices for 20 subjects, 3 temporal scales (encompassing 1 experiment, 3 experimental sessions, and 
25 time windows), and 1 random network null model. In the dynamic connectivity investigation, we 
constructed connectivity matrices for 20 subjects, 18-34 time windows, 3 different null models, 10 values 
of the inter-layer coupling C, and 4 values of time window length (80, 90, 100, and 110 TRs). In light 
of the computational extent of this work, we note that we did not employ Kernighan-Lin (KL) node- 
swapping steps [61] in our optimization of Q pm or Q m i, as they would be computationally prohibitive and 
are not necessary in the present context. KL steps move individual nodes between communities in order 
to further optimize a single sample of Q pm or Q m i [33,62,63]. As we focus on the mean properties of 
ensembles of partitions (and use them to report reliable measurements of architectural properties) rather 
than on the values of diagnostics for any individual partitions, KL steps that provide a marginal increase 
in the value of Q pm or Q m i would not be helpful for our study. 
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Precentral gyrus 


Parahippocampal gyrus, posterior 


Temporal pole 


Lingual gyrus 


Superior temporal gyrus, anterior 


Temporal fusiform cortex, anterior 


Superior temporal gyrus, posterior 


Temporal fusiform cortex, posterior 


Middle temporal gyrus, anterior 


Temporal occipital fusiform cortex 


Middle temporal gyrus, posterior 


Occipital fusiform gyrus 


Middle temporal gyrus, temporooccipital 


Frontal operculum cortex 


Inferior temporal gyrus, anterior 


Central opercular cortex 


Inferior temporal gyrus, posterior 


Parietal operculum cortex 


Inferior temporal gyrus, temporooccipital 


Planum polare 


Postcentral gyrus 


Heschl's gyrus 


Superior parietal lobule 


Planum temporale 


Supramarginal gyrus, anterior 


Supercalcarine cortex 


Supramarginal gyrus, posterior 


Occipital pole 


Angular gyrus 


Caudate 


Lateral occipital cortex, superior 


Put amen 


Lateral occipital cortex, inferior 


Globus pallidus 


Intracalcarine cortex 


Thalamus 


Frontal medial cortex 


Nucleus Accumbens 


Supplemental motor area 


Parahippocampal gyrus (superior to ROIs 34,35) 


Subcallosal cortex 


Hippocampus 


Paracingulate gyrus 


Brainstem 



Table 1: Brain regions present in the Harvard-Oxford Cortical and Subcortical Parcellation Scheme pro- 
vided by FSL [2,9]. 
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Type of Flexibility 


Session 


Cortical 


Connectional 


Nodal 


Temporal 


F 














1 

2 
3 


0.027±0.009 
0.030±0.009 
0.027±0.008 


0.041±0.008 
0.042±0.007 
0.039±0.008 


0.070±0.001 
0.067±0.001 
0.064±0.001 


0.027±0.010 
0.030±0.009 
0.026±0.009 


F' 














1 

2 
3 


1.93±0.30 
1.95±0.27 
1.86±0.25 


2.58±0.25 
2.53±0.23 
2.43±0.23 


3.19±0.01 
3.01±0.01 
2.92±0.01 


1.90±0.31 
1.94±0.28 
1.85±0.26 



Table 2: Flexibility in cortical network and random network null models for two different definitions of flexibility. 
One definition is F = 2fcLi fit wnere fi ls tne number of times node i changes module allegiance in a given experimental 
session divided by the total number of possible changes (i.e., by the number of consecutive pairs of layers in the multilayer 
framework). The second definition is F' = 2»=i fl> where is defined as the total number of modules of which node i 
is a member at some point in the experiment session. We also indicate mean values and standard deviations. 



Type of Flexibility 


Session 


Null Model 


^-statistic 


p- value 


F 


1 


Connectional 


-10.62 


1.2 x 10" 18 






Nodal 


-47.56 


1.1 x 10" 75 






Temporal 


-0.64 


0.51 




2 


Connectional 


-9.47 


5.5 x 10 -16 






Nodal 


-43.41 


1.7 x 10" 71 






Temporal 


-0.16 


0.86 




3 


Connectional 


-9.86 


1.1 x 10" 16 






Nodal 


-46.95 


4.5 x 10 -75 






Temporal 


0.77 


0.43 


F' 


1 


Connectional 


-14.54 


1.8 x 10" 27 






Nodal 


-43.44 


1.6 x 10 -71 






Temporal 


4.67 


8.4 x 10" 6 




2 


Connectional 


-14.04 


2.3 x 10" 26 






Nodal 


-39.67 


2.1 x 10" 67 






Temporal 


3.95 


1.3 x 10 -4 




3 


Connectional 


-14.87 


3.5 x 10" 28 






Nodal 


-43.80 


6.9 x 10" 72 






Temporal 


2.79 


0.0061 



Table 3: Results for one-sample t-tests on flexibility in cortical networks versus random network null models 

for two different definitions of flexibility (F and F'). We report these results for the connectional, nodal, and temporal 
null models. Negative t-statistics indicate that the flexibility of the null model network is greater than that of the cortical 
network, and positive t-statistics indicate that the flexibility of the cortical network is greater than that of the null model 
network. 
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Figure 1: Experimental Setup and Learning (A) Schematic of the cued sequence production (CSP) 
task. The response or "button" box (left) had four response buttons that were color-coded to match the 
notes on the "musical staff" (right) presented to the subject in the visual stimulus. This visual stimulus 
was composed of 12 notes in sequence. Here we show one example of a single sequence. (B) Movement 
time as a function of practiced trials, whose decreasing slope indicates that learning is occuring. (We 
have aggregated trials into 10 trial bins per session.) 
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Figure 2: Properties of the static and dynamic modularity indices Q± and Q m i- The mean 
(column 1), standard deviation (column 2), and maximum (column 3) of the static modularity index 
Q± is shown for (A) the large scale (entire experiment), (B) the mesoscopic scale (three experimental 
sessions), and (C) the small scale (individual time windows) over the 100 samplings. Row (D) shows the 
mean (column 1), standard deviation (column 2), and maximum (column 3) of the dynamic modularity 
index Q m i over the 100 samplings. In the figure, the standard deviation is abbreviated as STD. Boxplots 
indicate 95% confidence intervals over subjects. 
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Figure 3: Effects of the coupling parameter C on the four diagnostics in this study: modularity 
index Q m i, number of modules n, module size (i.e., number of nodes) s, and module stationarity £. We 
first averaged values over 100 'optimal' partitions (see the discussion in the text), so this figure gives 
mean values of all diagnostics. The error bars indicate standard deviations over subjects and sessions. 
Colors indicate network type: cortical network (blue), temporal null model network (green), nodal null 
model network (orange), and connectional null model network (red). Error bars for different network 
types at a given value of C (0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2) are offset from each other for better 
visualization. 
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Figure 4: Effect of the time window length T on the four diagnostics in this study: modularity index 
Qmi, number of modules n, module size (i.e., number of nodes) s, and module stationarity £. We first 
averaged values over 100 'optimal' partitions (see the discussion in the text), so this figure gives mean 
values of all diagnostics. The error bars indicate standard deviations over subjects and sessions. In the 
figure, we give time windows in terms of the number of data points in the time series (i.e., the number 
of TRs). Colors indicate network type: cortical network (blue), temporal null model network (green), 
nodal null model network (orange), and connectional null model network (red). Error bars for different 
network types at a given value of T (80, 90, 100, 110) are offset from each other for better visualization. 
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